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C/3 Purpose: Current inverse planning methods for IMRT are limited because they are not de- 

signed to explore the trade-offs between the competing objectives between the tumor and normal 
tissues. Our goal was to develop an efficient multiobjective optimization algorithm that was flex- 

£>~> ible enough to handle any form of objective function and that resulted in a set of Pareto optimal 

-G plans. 

Methods: We developed a hierarchical evolutionary multiobjective algorithm designed to quickly 
generate a small diverse Pareto optimal set of intensity modulated radiation therapy (IMRT) plans 
that meet all clinical constraints and reflect the optimal trade-offs in any radiation therapy plan. 

^ The top level of the hierarchical algorithm is a multiobjective evolutionary algorithm (MOEA). The 

t-H genes of the individuals generated in the MOEA are the parameters that define the penalty func- 

tion minimized during an accelerated deterministic IMRT optimization that represents the bottom 
level of the hierarchy. The MOEA incorporates clinical criteria to restrict the search space through 
protocol objectives and then uses Pareto optimality among the fitness objectives to select individu- 
als. The population size is not fixed, but a specialized niche effect, domination advantage, is used 
which controls the patient size and plan diversity. The number of fitness objectives is kept to a 
minimum for greater selective pressure, but the number of genes is expanded for flexibility that 
allows a better approximation of the Pareto front. 

Results: Acceleration techniques implemented on both levels of the hierarchical algorithm resulted 
• i-h in short, practical runtimes for multiobjective optimizations. The MOEA improvements were eval- 

uated for example prostate cases with one target and two OARs. The modified MOEA dominated 
11.3% ± 0.7% of plans using a standard genetic algorithm package. By implementing domination 
advantage and protocol objectives, small and diverse populations of clinically acceptable plans 
that were only dominated 0.2% ± 0.05% by the Pareto front could be generated in a fraction of 
an hour. 

Conclusions: Our MOEA produces a diverse Pareto optimal set of plans that meet all dosimetric 
protocol criteria in a feasible amount of time. Our final goal is to improve practical aspects of the 
algorithm and integrate it with a decision analysis or human interface for selection of the IMRT 
plan with the best possible balance of successful treatment of the target with low OAR dose and 
low risk of complication for any specific patient situation. 
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1 Introduction 



Intensity modulated radiation therapy (IMRT) uses inverse planning optimization to determine com- 
plex radiation beam configurations that will deliver dose distributions to the patient that are better 
than could be achieved using 3D conformal radiation therapy. The problem is inherently multiob- 
jective with competing clinical goals for a high, uniform target dose and low dose to the organs 
at risk (OAR). The IMRT inverse planning algorithms commonly optimize a single objective func- 
tion (also known as "penalty function") that is a linear combination of many separate organ and 
tumor-specific objective (penalty) functions. The solution of this problem is one point in the space 
of feasible treatment plans. To find another point in this space, the user can vary any one of the 
parameters that make up the penalty function. Current clinical algorithms do not provide methods 
for efficiently searching this space of feasible plans. The search is complicated by the fact that in 
comparing plans, one plan may have achieved superior values for some objectives but inferior values 
for other objectives. 

A multi objective optimization (MOO) algorithm is one method for searching through the space of 
feasible plans. Such an algorithm varies the parameters that define an objective function and cal- 
culates a new plan for each new function. For each plan created, the algorithm calculates objective 
values that characterize its quality. If all of the objective values of one plan are better or the same as 
the values of a second plan, the first plan is said to dominate the second. A plan is said to be Pareto 
optimal if there are no plans that dominate it. Note that any Pareto optimal plan has some superior 
objective values compared to other Pareto optimal plans, but inferior values for others. MOO algo- 
rithms must generate a large number of plans that each require an individual IMRT optimization. 
The resulting computation time could make this algorithm impractical depending on the efficiency 
of IMRT optimization method used and the number of plans required. The final product of a MOO 
algorithm is a set of Pareto optimal plans, known as the Pareto front, so ultimately, there must be a 
method to select the individual plan to be used. 

The true Pareto front of this multiobjective problem is a multidimensional surface that contains an 
infinite number of points, and hence there is the potential for the production of Pareto populations 
that are too large for practical purposes. Methods for handling this issue, called niche functions, 
come in many forms, but are all designed to limit the number of individuals that can survive in a 
population while maintaining diversity. For stochastic multiobjective algorithms to be practical, the 
use of an appropriate niche function may be quite important. 

The nature of the separate objective functions used in the MOO plays an important role in the 
optimization. Those functions that are convex may be optimized using deterministic algorithms. 
Stochastic optimization methods are needed for optimizing those functions that are not convex to 
avoid becoming trapped in local minima. In selecting a multiobjective optimization approach, the 
user may have to make a decision between using a deterministic method that is fast but limits the 
type of objectives and a stochastic method that is very flexible but slower. If different objective 
functions are used to evalute the quality of dose distributions than are used in the single penalty 
function that optimizes beamlet intensities, it is possible that some parts of the optimization can be 
done deterministically while other parts may use stochastic methods. 

One type of stochastic optimization algorithm is an evolutionary algorithm which uses concepts 
based on genetics and natural selection to create a population of new solutions and to select the 
fittest among them. An "individual" is described by a set of genes. New individuals are created 
through recombination of genes from two "parent" individuals selected from the existing population. 
Random mutations can occur which help prevent the population from becoming trapped in local 
minima. Natural selection of superior individuals followed by reproduction using their genes leads 
the algorithm to find the Pareto optimal solutions. In practice, further heuristic methods are often 
needed to speed up this process. 
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This paper describes the development of a multi objective optimization algorithm that combines 
an evolutionary method to search feasible plan space with a deterministic method using a single 
objective penalty function that quickly finds the optimal beamlet intensities for that function. While 
the penalty function must be convex, the hierarchical method allows us to use any fitness objective 
(convex or non-convex) to evalute plan quality without losing too much in speed. Utilization of a 
MOO algorithm requires the selection of one plan from the Pareto optimal set to be used for IMRT 
treatment of the patient. There are many methods for doing so. We have selected an a posteriori 
method of multiobjective optimization where plan selection is done after the Pareto front is found 
by either an influence diagram IT] O or inspection by a radiation oncologist. With both of these 
approaches, it is best to restrict the set of plans to a relatively small number, e.g. 10-20, and this is 
reflected in our MOO algorithm. 

In the methods section, the hierarchical algorithm is described with the emphasis on the evolutionary 
element and how it interacts with the deterministic step. In the results, performance of the algorithm 
is characterized using two clinical prostate cancer cases. 

2 Methods and Materials 

A multiobjective evolutionary algorithm (MOEA) has been developed for the top level of a hierar- 
chical IMRT optimization method. The lower level is a deterministic algorithm that incorporates a 
constrained quadratic optimization algorithm (BOXCQP). We have used a hierarchical process be- 
cause there are two sets of variables we would like to optimize: (1) beamlet intensities (determin- 
istic optimization), and (2) the objective parameters such as maximum dose, dose homegeneity etc 
(stochastic optimization) . At the stochastic level of the algorithm, an individual is created and each 
individual is characterized by a set of genes that determine the relative weighting between OARs 
and targets and a set of genes that describe the OAR penalty functions. The deterministic algorithm 
incorporates these genes into a quadratic objective function which is then solved for the beamlet 
intensities. The MOEA uses a separate set of objectives to determine the quality of each individual 
treatment plan. As mentioned above, considerations of speed have led us to incorporate several 
heuristic stages in our algorithm which results in a faster convergence to a clinically-acceptable set 
of Pareto optimal plans. 

The following sections provide more detail of the algorithm as well as descriptions of the methods 
used to characterize the performance of the MOEA. 

2.1 Multiobjective Evolutionary Algorithm 

An evolutionary algorithm was chosen as the stochastic component of the multiobjective method 
because of its effectiveness in searching large regions of solution space while still incorporating 
techniques for focusing the search on regions more likely to produce better plans 131 |4j [5j |6l [7J . A 
number of characteristics of the IMRT optimization problem guided the development of our MOEA. 
The critical factors were (1) the very large search space, (2) the relatively large number of objectives, 
(3) the need for solutions to conform to "clinically-acceptable" criteria, and (4) the time constraints 
for clinical use. In short, our modified MOEA was designed for faster and more effective convergence 
to the Pareto front inside the clinically acceptable space. The elements of the MOEA and highlights 
of the impact of these critical factors are described below. 

2.1.1 Enhancements to the MOEA 

Our MOEA uses no fixed number of individuals per generation, no summary single fitness function, 
and no tournament style selection method (as are common in traditional genetic algorithms). In- 
stead, individuals are created one at a time and are eliminated if they are Pareto dominated by any 
existing individual or if they violate any clinical constraints. Conversely, if the new individual domi- 
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nates any of the current population, those dominated individuals will be eliminated. If the individual 
survives, its genes become candidates for recombination in the next individual. In addition to these 
algorithmic modifications, various parameters of the multiobjective method, such as mutation rate, 
have been optimized for performance. 

We also implemented a specialized niche function, called domination advantage. When calculating 
whether plan A dominates plan B, each objective value of plan A only needs to be within an e of 
the corresponding plan B objective in order to dominate it. This domination advantage, e, can 
be thought of as the difference in objective values below which the clinical significance cannot be 
determined. Thus, it provided a means to avoid clustering of solutions so that the Pareto optimal set 
spans the entire clinically relevant space. It was also used to limit the total number of plans. The 
e used in domination advantage was not a constant but was set equal to c • (n current - n goal ). The 
constant c was adjusted so that the final number of plans was approximately the desired number, 

^current ^ ^goal- 

If two individuals were being compared to test for domination and were very close together, i.e. 
with all objective values within e of each other, both individuals would 'dominate' the other using 
domination advantage. In this rare case, a single fitness function was used that was a weighted sum 
of fitness objectives using importance weighting predetermined by the user. The individual with the 
worse single fitness function would be eliminated and the other would survive. Note that the result 
is different from might occur in traditional evolutionary algorithms where two individuals from 
different regions of the objective space would be forced into such a comparison. In our algorithm, 
these two individuals are so close together, the assumption is that little information of clinical value 
would be lost by the undesirable, but necessary use of a single fitness function elimination. 

We used two sets of objectives in the stochastic part of the algorithm, which we call fitness objectives 
and protocol objectives. Fitness objectives are the dosimetric functions that best represent the clinical 
goals and are used to determined plan quality. Protocol objectives are functions that delineate certain 
dosimetric cutpoints that need to be met, usually because of clinical protocol restrictions. Protocol 
objectives will restrict the search space such that the final set of plans are guaranteed to meet all 
clinical criteria. As will be described in more detail below, these two sets of objective functions are 
distinct and play different roles during the optimization. 

2.1.2 Fitness and Protocol Objectives 

In multiobjective optimization, smaller numbers of fitness objectives tend to result in greater selec- 
tive pressure and better performance. Thus, it is desirable to have only one fitness objective for 
each structure. We chose effective uniform dose (EUD) as an appropriate objective for OARs; other 
functions such as complication probability functions could also be used. Unlike for OARs, we did 
not find target EUD to be an effective objective, because a uniform target at prescription dose could 
be deemed equivalent to a dose distribution with high variance and a higher mean dose. Typically, 
target dose distributions must meet minimum dose requirements as well as limits on inhomogeneity. 
We removed the minimum dose objective by prioritizing it and scaling the dose distributions to meet 
all minimum dose requirements, thereby leaving maximum dose (dose homogeneity) as the sole 
target objective. 

The two sets of objective functions that were used are: 
• Fitness objectives 
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ff%t{x) = (//it,l(x), • • • , ffit, m [x), ffit,m+l(x), ■ ■ ■ , ffit,m+n{%)), 

where m = # of target fitness objectives 

n = # of OAR fitness objectives 
ffit,%{%) = d? ax for i e {1, . . . , m}, i.e. for targets 

nvoxi 

fm,%{x)= EUD(d) = (l/nvoxi df) 1/a "- for i e {m + 1, . . . , m + n}, i.e. for OARs. 

(1) 

where is the dose to j th OAR voxel, nvoxi is the number of voxels in the i th OAR, c?™ ax is the 
maximum dose to any voxel within the ith target, and <n is a constant (> 1) that characterizes 
the radiation response of each of the n OARs J5J |9j HUH . 

Protocol objectives 

fprot{x] — \fprotfX\Xji' ' ' ) fprot,p\&)t fprot,p+l\%)) • • • j fprot,p+q\^))i 

where p = # of target protocol objectives 
q = # of OAR protocol objectives 

k k pro *' for all clinical maximum dose contraints 

otherwise 

{Tivox if nvoXviolate > k ^> y 

mo a e ' nvox k £ or a jj c li n ±cal dose volume constraints 

otherwise 

(2) 

where D prot ^ is the maximum dose allowed allow to the structure associated with the k th 
protocol objective, nvox V i i ate ,k is the number of voxels in the OAR associated with the k th 
protocol objective that receive doses that exceed the allowed dose of a dose volume constraint 
(DVC), nvox k is the total number of voxels in the OAR associated with the k th protocol objec- 
tive, and VoAR,k is the maximum fraction of the OAR associated with the k th protocol objective 
allowed over the dose for a DVC. Minimum target dose constraints are achieved by scaling the 
dose distribution if necessary. 



2.1.3 Stages of the MOEA 

There are three stages of the MOEA that correspond to: (1) a broad search of the feasible plan 
space, (2) a focusing of the search to a region of space that meets protocol requirements and (3) 
a search of the feasible space within the protocol constraints. Throughout the last two stages, the 
domination advantage procedure is used to both spread the plans out (avoiding clustering near local 
minima) and to limit the number of plans. As the optimization proceeds, the current set of solutions 
approaches the true Pareto front defined by the fitness objectives. 

In stage 1, the individuals were created by random selection of genes. The genes were used to 



construct the quadratic objectives of the deterministic part of the algorithm (see Sec. 2.2). Then 
the fitness and protocol objective values were calculated, and the individual was added to or re- 
jected from the current population, whose size was n current . When n current equaled the pre-defined 
number of desired plans determined by the decision procedure, n goa i, stage 1 ended. 

In stage 2, new individuals were created by recombination of current genes and mutation. During 
this stage, the direction of mutations were guided depending on whether or not the parent plan met 
the protocol objectives in order to more quickly search relevant space. Selection proceeded as in 
stage 1 with the addition of the domination advantage procedure. In this stage, the population was 
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allowed to grow larger than n goa i. Stage 2 ended when there were n goa i individuals who met all of 
the clinical criteria (protocol objective values = 0), at which time, all individuals who did not meet 
these objectives were eliminated. 

Stage 3 proceeded similarly as stage 2 except that the protocol objectives were converted to con- 
straints, mutation direction was not guided, and Pareto optimality was determined only by the fitness 
objectives. This stage ran for a pre-determined number of iterations. 

One reason to separate the objectives into two sets is to reduce the dimensionality of the final 
solution space. As the number of fitness objectives considered in Pareto domination increases, the 
probability of domination is reduced by a factor on the order of ~ .5«o6i«cti»«i j an d hence selective 
pressure was reduced. By converting protocol objectives into constraints, the solution space was 
restricted to a smaller number of dimensions equal to the smaller number of fitness objectives. 

A multi objective optimization algorithm needs to be paired with a decision making component. This 
will have an impact on the number and variation of the solution set. In our case, the decision maker 
is either a human or an influence diagram, neither of which can handle large numbers of solutions 
or very nearly identical solutions. Hence, our algorithm was designed to produce on the order of 
n goa i plans, i.e. 10-15 plans for n goa [ = 10. The algorithm can be adjusted to produce any size 
solution set. 



2.2 Deterministic Algorithm 

For the beamlet intensity optimization, we employed a modification of the Bound Constrained Con- 
vex Quadratic Problem (BOXCQP) algorithm [TTT1 fT2l [TBI . This efficient method was implemented 
to minimumize a quadratic objective function with the bound constraint that all beamlet intensities 
are non- negative and less then some maximum allowed beamlet intensity, x max . Non-zero penalty 
factors were assigned to OAR voxels that had exceeded the OAR dose parameter in each iteration. 
A short description is provided. 

To compute optimal beamlet intensities, BOXCQP minimized the matrix form of the penalty function 
defined by each individual's genes summed with a smoothing function as described below: 

minimize f det {x) = (d- D det ) T P(d- D d et)\j =B . s 

+ k(Sx) t (Sx) (3) 
subject to < x b < x max 



where f det is the objective function for the determinstic component of the algorithm, x is the vector 
of the beamlet intensities, B is the precalculated matrix storing the dose to each voxel from each 
beamlet, P is a diagonal matrix of all importance weighted penalty factors for each voxel, k is the 
smoothing constant, S is the smoothing penalty matrix, and Xb represents any element in vector x 

A smoothing term was introduced because BOXCQP depends on solving a well-posed problem Hill 
[T2ll . It is also desirable to find the deliverable solution for a beamlet weight vector x that is not 
excessively modulated due to delivery limitations, i.e. MLC movement. The smoothing term was 
based on the fact that the smaller the second derivative of a fluence map, the easier it is to deliver 
in practice. Therefore, S was the second derivative matrix by finite difference method, and k was a 
smoothing coefficient that was found empirically. 

The algorithm was implemented using MatlatQ Since this code was executed many times during 

the course of the optimization, consider able effort was taken to speed execution. 

ifhe MathWorks, Nftick, MA r 
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2.3 Treatment planning platform 



A treatment plan consists of a set of beams, which describe the geometry of radiation delivery. Each 
radiation beam was divided into a 2D array of beamlets. Each beamlet represents a .6 cm by 1 cm 
area of the beam. A beamlet length of 1 cm was selected because of the size of the collimator leaf, 
and a width of 0.6 cm was used as a compromise between speed and resolution. Beamlet intensities 
were the variables optimized in the deterministic IMRT optimization algorithm. 

The patient was considered as a collection of 3D voxels that are segmented into a set of targets and 
critical normal tissues (OARs). A model of radiation transport mapped each beamlet intensity to the 
dose it deposited at each voxel; this was stored in a beamlet-to-voxel dose contribution matrix, which 
was calculated once before the optimization. Typically the number of beamlets ranged from 1000- 
5000 while the number of voxels ranged from 10 4 to 10 5 . The entire system was modeled using 
the Prism treatment planning system H141 PT5l [TBI that was augmented with a pencil beam dose 
calculation algorithm H17H . The optimization algorithms interacted with Prism via file input/ output. 

2.4 Experiments 

As the focus of this work is to introduce this multi objective method, results are shown for only two 
prostate cancer example cases. Case 1, used in all experiments, was selected to be challenging and 
had a PTV that was ~8.5 cm across at the widest point and had 39% overlap with the bladder and 
32% overlap with the rectum. Case 2, used only in sections 3.5 and 3.6, had a PTV that was ~7cm 
in diameter with only 9% overlap with the bladder and 17% overlap with the rectum. All plans used 
seven equally spaced 6 MV photon beams. 

Three fitness objectives were used in all of the optimizations: bladder (or bladder wall) EUD, rectal 
wall EUD, and target maximum dose. OAR voxels inside the PTV were excluded from objective 
calculation. All dose distributions were scaled to meet minimum target dose requirements, thereby 
making the target maximum dose a measure of target variance. The protocol objectives represented 
a maximum allowed target dose and OAR DVH constraints. A comparison was also made for two 
example prostate patients between the MOEA and a plan generated using a commercial objective 
deterministic IMRT optimization algorithm with the current clinical protocol. Of all possible MOEA 
plans, we chose the one plan from the MOEA output that most closely matched the clinical plan's 
target dose distribution. 

It should be noted that stochastic algorithms do not guarantee convergence to the Pareto optimal 
set. The degree to which this algorithm approximates that set is the subject of several of the result 
sections described below. 

The following experiments were performed to evaluate individual aspects of the algorithm: 

• demonstration of the evolution of a population of IMRT plans towards the Pareto front through 
the three stages of a multiobjective optimization (Section 3.1); 

• a mulitobjective comparison of performance between a standard genetic algorithm package 
and the modified evolutionary algorithm (Section 3.2); 

• demonstration of the role of domination advantage in controlling population size, maintaining 
diversity, and improving performance (Section 3.3); 

• evaluation of the performance of the MOEA with and without protocol objectives (Section 
3.4); 

• a DVH comparison between plans produced by the MOEA and plans produced using our cur- 
rent clinical protocol and a commercial planning system (Section 3.5); 

• assessment of the ability of the MOEA to approximate the true Pareto front (Section 3.6). 
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A common method for comparing two different algorithms is to select the 'best' plan from each 
population and do DVH comparisons; however, such a method is not a multiobjective metric that 
can compare algorithms that produce populations of plans. It focuses on only a small subset of the 
results, does not lead itself to quantitative or statistical analysis, and, given the multiobjective nature 
of the problem, selection of the plan used for comparison is subjective. As previously discussed, in 
multiobjective optimization a plan can only be deemed superior to another if it is better in all of its 
objectives, which is the definition of Pareto domination. In our evaluation of different multiobjective 
algorithms, we used a metric based on Pareto dominance that we term "domination comparison", 
Cd- This is the percentage of all possible comparisons between two groups of plans in which the 
plan from the first population dominates the plan from the second population. In mathematical 
form, the domination comparison metric for a population resulting from one algorithm, U, in com- 
parison with a population of plans obtained with another algorithm, V, is defined as follows: 

Let the set of plans obtained by formulation U = u 
and those obtained by formulation V = v 

and Cd = \u\*\v\ (4) 

where f(ut, vj ) = 1 if dominates Vj 
and f(ui,Vj) = otherwise. 



Domination comparison is a purely multiobjective metric that gives some indication of which of two 
populations is closer to the Pareto front; however, it does not give a measure of the diversity of 
the population or coverage of the Pareto front. For this reason we included 3-D Pareto plots for 
graphical display of the final populations for the two prostate examples. For greater ease, we made 
all of these plots uniform, with rectum objective on the x-axis, bladder objective on the y-axis, and 
target objective in color scale. 



3 Results 

3 . 1 Stages of the MOEA 

A population of plans for case 1 was tracked through the three stages of the evolutionary algorithm 
(Figures la, lb, and lc). Three fitness objectives were used: rectal wall EUD, bladder EUD, and 
maximum target dose (MT). Only OAR voxels outside the target were used in the calculation of the 
fitness objectives. 

Often Pareto optimality is demonstrated using only two objectives which makes determining the 
Pareto front visually easy. Given that we used three fitness objectives, we used color to indicate the 
value of the third axis. For uniformity, rectal wall EUD was plotted on the X-axis, bladder EUD on 
the Y-axs, and the value of the maximum target dose (MT) was represented in color scale. The plots 
illustrate the difficulty in visualizing even a 3-dimensional fitness objective space, which is small 
compared the the dimensions of most IMRT planning. 

Dose distributions were scaled to meet all minimum target dose requirements to eliminate the need 
for a minimum target dose fitness objective. The maximum target dose represents clinically unde- 
sirable hot spots and is a measure of target variance. Throughout the presentations of these results, 
lower fitness objective values correspond to more optimal plans. We used n goa i = 10 throughout 
this work. 

The population after the first 10 iterations obtained from completely random genes have relatively 
poor objective values. Stage 2 is complete after 53 iterations, and objectives values have improved. 
The ten clinically acceptable plans are shown with a black box outline. All other plans were removed 
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before stage 3 began. For this example, the greatest amount of improvement is seen in a small 
number of iterations during stage 2 as the initial population was not close to Pareto optimal. Stage 
3 was completed after a total of 200 iterations. Further improvement is observed in the objective 
space during stage 3. 

3.2 Performance of the Hierarchical Evolutionary Algorithm 

The modified evolutionary algorithm was compared to the Matlab standard genetic algorithm pack- 
age that used a single stage, a set population size, and used a tournament- style selection method 
with a single fitness function. Both algorithms used the previously described deterministic algo- 
rithm. To isolate differences in the algorithms, domination advantage and protocol objectives were 
not used in either multiobjective optimization. Important differences included the selection method, 
allowed population size, and mutation rates. Bladder EUD, rectal wall EUD and maximum target 
dose (MT) were used as the fitness objectives. 

Twenty-five multiobjective optimizations that each produced a set of 11 to 15 plans were performed 
for both situations to establish statistical significance, and each multiobjective optimization executed 
a total of 200 BOXCQP iterations in ^10 minutes. Objective values (OAR EUD and target maximum 
dose) for individual plans resulting from two optimizations are shown in Figure 2. 

Populations resulting from our MOEA dominated populations resulting from the standard genetic 
algorithm D c = 11.3% ± 0.7% of possible comparisons versus D c = 0.04% ± 0.02% in reverse. 
Results indicate that the new evolutionary algorithm was able to produce a population of plans with 
lower OAR EUD values and lower target variance in a short amount of time. 

3.3 Domination Advantage 

Our MOEA does not have a set population size or limit and does not use a single-objective fitness 
function for selection. Resulting populations could potentially be very large even with a relatively 
small number of objectives. To keep the population size in the desired range, to increase selective 
pressure, and to promote diversity in the population, a specialized niche effect, domination advan- 
tage, was used. With domination advantage, once the population was greater than n goa i, potentially 
dominating individuals were given an advantage when testing for domination. One or more objec- 
tives of a dominating individual can be inferior as long as it was within e of the objective of the 
individual being tested. The size of e increased linearly with the difference between the n current 
and n goa i. A comparison of two multiobjective optimizations of a prostate case (with and without 
domination advantage) was performed to demonstrate the effect on a single population (Figure 3). 
Twenty-five multiobjective optimizations were performed with and without domination advantage 
to establish statistical significance. 

Average population size without domination advantage was 68.7 ± 2.1 versus 13.2 ± 0.2 with 
domination advantage. The population generated using the niche effect dominated an average 
of D c = 1.73% ± 0.35% of the population generated without it and D c = 0.17% ± 0.04% of 
dominations occurred in reverse. In addition to creating a smaller, diverse population, improvement 
in the resulting population was observed when domination advantage was used. When domination 
advantage was used, there were fewer individuals for a more managable final population size, and 
individuals were closer to the true Pareto front possibly due to the increase in selective pressure from 
the niche effect. Domination advantage did not allow searching in regions that give only a small 
improvement in one objective at great cost to other objectives. This speeds up the optimization; 
however, it also means that some of the Pareto front may not be mapped. 
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3.4 Effect of Protocol Objectives 



Protocol objectives were used to focus the search to a region of objective space that met pre- 
determined clinical dosimetric requirements. Performance of multiobjective optimizations were eva- 
luted for an example prostate case with and without protocol objectives that reflected constraints 
from a clinical protocol. The effectiveness of the preference space guidance algorithm is shown in 
Figure 4. 

In terms of Pareto domination the protocol objectives did not improve or degrade results; however, 
without protocol objectives, the search focused on too wide a search space. Of all individual plans 
generated without using protocol objectives, only 4.4% ± 0.8% met the clinical criteria, and only 
68% of the full multiobjective optimizations had at least one plan that met all clinical requirements. 
In other words, 32% of the optimization procedures did not yield a usable plan when the protocol 
constraints were considered, whereas all of the plans generated using the protocol objectives met all 
of the clinical criteria. 

3.5 Two Clinical DVH Comparisons 

The MOEA was compared to our clinical planning methods for the two prostate cases. IMRT plans 
were optimized using a 200 iteration MOEA optimization and using our standard clinical inverse 
planning procedurq^J The MOEA plan whose target dose distribution most closely matched the plan 
using the clinic's deterministic algorithm was chosen from each multiobjective population for each 
case. A DVH comparison was made between these plans and the corresponding clinical plans. None 
of the example patient studies included femur and unspecified tissue objectives, and the MOEA plans 
were not translated into deliverable beams. This comparison is only a demonstration of the potential 
of this algorithm and not designed to be a comparison of plans ready to be used in the clinic. 

The DVH comparisons for the target, rectal wall, and bladder for case 1 are shown in figures 5a-5c. 
The minimum target dose requirements were matched for these plans; however, the uniformity of 
the target dose distribtuion and lower maximum target dose were better for the MOEA plan. As this 
patient was more difficult especially with regards to the proximity of the entire bladder with the 
target, improvement of the MOEA OAR distributions over the clinical plan was modest. 32% of the 
volume of the rectum was inside the target, and some of this volume received higher dose in the 
MOEA plan. Once outside the PTV, the dose to the rectum dropped off more quickly for the MOEA 
than for the plan produced using the current clinical protocol. 

The DVH comparisons for the target, rectal wall, and bladder wall for case 2 are shown in figures 5d- 
5f. For this prostate patient example, there was more separation between the OARs and the target, 
and this resulted in a greater potential for improvement using the MOEA. The target distributions 
nearly overlap, with the clinical plan's coldest 1% of voxels below 75.5Gy versus 75.4Gy for the 
MOEA plan, and the clinical plan's hottest 1% of voxels above 79.1 Gy versus 79.3 Gy for the MOEA 
plan. Dose levels for both the bladder and especially the rectum are much lower for the MOEA plan 
over the entire DVH. 

3.6 Estimation of the Pareto front 

The same two cases were used to investigate the ability of the MOEA to approximate the true Pareto 
front. For case 1, a long multiobjective optimization generated 9118 individuals in 12.6 hours 
without domination advantage in order to approximate the true Pareto front. Results were compared 
to 25 shorter optimizations that each generated 200 individuals using domination advantage. Figure 
6a plots the results from one of the 25 short runs performed using domination advantage against 
the long optimization without using it. A similar comparison using case 2 was also performed using 

2 Pinnacle 3 , Philips Medical Systems 
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both a long (5000 individuals) run and 25 short runs (200 individuals) each. Results from one of 
the short runs is plotted with the long run results in Figure 6b. 

For case 1, each individual in the long run dominated only Dc = 0.20% ± 0.05% of the popula- 
tion generated using the short optimization. For case 2, each individual in the true Pareto front 
dominated an average of only Dc = 0.14% ± 0.04% of the population generated using the short 
optimization. This demonstrates how quickly the algorithm with domination advantage approaches 
the Pareto front. The smaller populations were diverse and covered a wide range of the Pareto front, 
but certain regions characterized by diminishing returns were missed when using the domination 
advantage feature. 

In order to see if the small performance improvement using domination advantage in short runs 
(Section 3.3) translated to longer runs, the long optimization runs described in the previous para- 
graph were repeated with and without domination advantage. For case 1, each individual in the 
long optimization using domination advantage dominated only Dc = 0.009% of the population 
generated using the long optimization without domination advantage with Dc = 0.004% in re- 
verse. For case 2 each individual produced using domination advantage dominated only Dc = 
0.005% of the population generated using the long optimization without domination advantage, 
and the comparison was Dc = 0.028% in reverse. Results suggest that for very long runs that 
closely approximate the Pareto front, domination advantage may slow convergence slightly, but dif- 
ferences are very small. Note that without the single fitness function domination advantage feature 
described in section 2.1.2, the MOEA would never get within ~ e of the Pareto front using domina- 
tion advantage, because all new individuals only slightly better than the existing population would 
have been eliminated. 

4 Discussion 

We have reported on the development of a hierarchical multi objective optimization algorithm that 
couples a stochastic evolutionary algorithm with a fast deterministic one. The multiobjective nature 
of the algorithm reflects our view that IMRT planning requires exploring possible tradeoffs between 
competing clinical goals. The decision making strategy that must be coupled with the the mul- 
tiobjective optimization has a profound influence on the structure of the algorithm. The current 
standard clinical decision making strategy supported by commercial inverse planning systems takes 
an interactive approach in which a plan is generated, evaluated by a human user, and a new search 
is made if necessary Unfortunately, these systems do not provide any method for performing an 
efficient and effective search of all possible tradeoffs between objectives. While the result is usually 
better than achievable with 3D conformal techniques, the decision maker is not aware of the large 
number of other feasible solutions that may be better. 

Another multiobjective method is the a posteriori approach in which a set of Pareto optimal plans 
are found first, and then this set is evaluated by a decision maker to select one. This is the nature of 
our MOEA and is a fairly popular approach for multiobjective problems in general H51 I18I . It has the 
advantage that the user does not have to determine relative importance amongst objectives before 
achievable tradeoffs are known. However, this comes at the price of time since there is little guidance 
in the search procedure. We have used several heuristic methods to focus the search. In addition, 
the domination advantage niche effect is designed to limit the final number of plans to the number 
that our decision process can handle. Bortfeld, Craft and colleagues have devised an alternative 
approach in which only convex objectives are used, the Pareto front is systematically mapped and 
calculation times are reduced by interpolating solutions Jl9j l20l 1211 1221 l23l . Specially designed 
software allows the decision maker to steer between solutions on the Pareto front and observe the 
trade-offs of the different regions. Another systematic search approach was taken by Lahanas and 
colleagues H241I25II in which attention was focused on the method of completely mapping the Pareto 
front with a simple decision making approach used a posteriori. An early paper by Yu which focused 
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on brachytherapy and radiosurgical planning presented an excellent discussion of multiobjective 
optimization in radiation therapy H26L He devised a plan ranking system that required some a 
priori importance input but used a simulated-annealing method to introduce some fuzziness and to 
broaden the search space. He also included a "satisficing" condition that functioned similarly to our 
protocol objectives. 

An a priori approach to multiobjective decision making in multiobjective optimization is the third 
alternative. Prior methods have the advantages of being among the fastest and providing a single 
solution. Jee et al [27 \ describe a method using the technique of lexigraphic ordering in which prior 
ranking of objectives leads to a sequence in which the first objective is satisfied and then turned into 
a constraint as the second objective is optimized and so on for the remainder. This is somewhat 
similar to our shifting the use of our protocol objectives between the second and third stages, but we 
make no preference between the objectives and all are satisfied simultaneously before transforming 
them to constraints. 

Our use of an evolutionary algorithm was determined by its extensive use in multiobjective opti- 
mization in other fields [[5j and because this type of algorithm can handle a much wider range 
of objective types than can deterministic algorithms. As with other aspects of our approach, the 
trade-off hinges on the speed of convergence. We used a hierarchical approach, coupled with the 
multiple stages, in order to mitigate this issue, and our results demonstrate that the performance is 
acceptable in a clinical setting. Wu and Zhu used a genetic algorithm for 3D conformal optimiza- 
tion to optimize the weighting factors J4|- Although they did not explicitly use Pareto optimality 
their choice of a genetic algorithm was based on a recognition of the multiobjective nature of the 
problem and a desire to be flexible with respect to single objectives. Evolutionary algorithms have 
found their greatest use for the solution of explicitly non-convex objectives, such as the problem of 
number and orientation of the radiation beams |f7 1l25[|28l|29H and the integration of leaf sequencing 
and intensity optimization H30II3U . Many other objectives in IMRT planning are also non-convex 
such as dose-volume objectives. For certain situations, it has been found that becoming trapped in 
local minima may be unlikely H321 or mathematical techniques are able to be applied to linearize 
non-convex objectives or approximate them with a convex surrogate H201I33L 

A critical component of any multiobjective optimization is the means by which two objects with 
multiple characteristics are compared or rated. A common example is the comparison of a dose- 
volume histogram from two different plans. When the DVHs do not cross, one can determine which 
is clearly better; when they do cross, "better" depends on the decision variables. This is the essence 
of Pareto optimality. In our IMRT optimization, we removed plans that were dominated within some 
epsilon since those plans were clinically similar or worse in every aspect when compared to other 
plans. When comparing algorithms or optimization methods, a similar multiobjective approach can 
be taken. In our evaluation of different versions of our algorithm, we found it more useful to base our 
definition of "better" on a statistical comparison using Pareto dominance as the defining criterion. 
A flexible multiobjective algorithm as described in this paper is the ideal means for making such 
comparisons which are more convincing than alternatives that are commonly used. 

5 Conclusion 

A multiobjective evolutionary algorithm was developed to find a diverse set of Pareto optimal and 
clinically acceptable IMRT plans. The novel selection method, flexible population size, and opti- 
mized mutation rates improved performance. The population diversity and size was effectively con- 
trolled using a niche effect, "domination advantage". The protocol objectives and multiple stages 
used in reproduction and selection guided the optimization to those regions of plan space that met 
clinical criteria. Key features of the algorithm are the ability to use any functional form of individual 
objectives and the ability to tailor the number and diversity of the output to the decision making 
environment. 
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Figure 1: Plots of the objective values of the plans at the end of (a) stage 1, (b) stage 2, and (c) 
stage 3. The colors represent different values for the target objective, MT. In (b), the square icons 
represent plans that met the clinical protocol objectives. 
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Figure 2: Plot comparing the objectives of the final populations of plans resulting from the hierar- 
chical approach with that of a standard genetic algorithm. The standard genetic algorithm is shown 
as asterisks, and the hierarchical approach is shown as solid circles. The color reflects the target 
maximum dose (MT), the target variance objective. Lower values in all objectives are more optimal. 
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Figure 3: Plot demonstrating the effect domination advantage has on a population of plans. The 
plans resulting from the hierarchical approach without using domination advantage are shown as 
asterisks, and those using domination advantage are shown as open squares. Maximum target dose 
(MT), the target variance objective, is displayed in color. Lower values in all objectives are more 
optimal. 
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Figure 4: Plot comparing multiobjective optimization with (WP) and without (NP) using protocol 
objectives. Bladder EUD and Rectum EUD are plotted for the sum of 25 full optimizations with 
and without using protocol objectives to restrict the search space. Maximum target dose (MT), the 
target variance objective, is plotted in color. The solid circles show the full Pareto front, and the 
open squares display the optimizations using protocol objectives. 
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Figure 5: DVH comparisons for bladder, bladder wall, rectal wall, and PTV are made between plans 
selected from the population produced by the MOEA and plans produced using the current clinical 
inverse planning method. 
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Figure 6: Plots for short multiobjective optimizations mapped onto an approximation of the Pareto 
front. Rectum and bladder EUD values are plotted with the target variance objective (MT) shown in 
color. The plans from the long optimization runs are represented as solid circles that make up most 
of the graphs. The plans resulting from the quick optimization are shown with asterisks using the 
same colorscale and surrounded by a black square for easier viewing. 
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